Distributed delays in systems of coupled phase oscillators 
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We study systems of identical phase oscillators with a delay distribution in their coupling that 
weights contributions arising from different past times. For any coupling topology with equal number 
of neighbors for each oscillator, we show that the frequency and stability of the fully synchronized 
states only depend on the mean of the delay distribution. However, transient dynamics leading to 
complete synchronization are affected by the shape of the delay distribution. Our work validates 
the use of discrete delays to study the synchronized states of more general delayed coupled systems. 
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Systems of coupled oscillators with delay in the cou- 
pling are important for several physical, chemical, engi- 
neering, and biological phenomena pQ. For example, cel- 
lular systems as the vertebrate segmentation clock [2]-[4] 
have been studied using systems of phase oscillators cou- 
pled with a discrete delay However, fluctuations 
in gene expression and transport of macromolecules in- 
troduce a variability in the time delay associated with 
cell to cell communication. One way to account for such 
variability is to consider that past states from different 
times can affect the evolution of the system with different 
weights. This is what is usually referred to as distributed 
delay 0. 

In this work we study a generic model of phase oscil- 
lators coupled with distributed delay: 
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where N is the number of oscillators in the system, 
k = 1, 2, . . . , N, 9k(t) is the phase of the fc-th oscillator, 
n is the number of coupling connections each oscillator 
has, K is the coupling strength, h(9) is a 27r-periodic 
coupling function, g(s) denotes the delay distribution, 
and ui is the intrinsic frequency of individual oscillators. 
The dki are the coefficients of the connectivity matrix D: 
dki = 1 if oscillator k is connected to oscillator I and 
dki = otherwise. We consider normalized delay distri- 
butions fulfilling J Q dsg(s) = 1, and define the mean of 
this distribution as: 
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where the phases of all oscillators are equal and grow lin- 
early in time with the global frequency fl. We substitute 
this ansatz into Eq. and obtain the solution for Q by 
a self-consistent approach [8j |9] : 
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Eq. Q is independent from the number of oscillators in 
the system. It depends only on the mean delay, and not 
on the particular shape of the delay distribution g(s). 
In particular, this means that the functional dependence 
of the global frequency in steady state is identical for 
discrete delays r and distributed delays with mean f — t. 

Linear stability of the synchronized steady state. Earl 
and Strogatz p]T] derived a stability criterion for the 
phase-locked steady state in systems of coupled phase os- 
cillators with discrete delays in the coupling, correspond- 
ing to g(s) — S(s — f) in Eq. 0. Their result holds for 
arbitrary coupling topologies where each oscillator has 
the same number of coupling connections, see Fig. [I] We 
extend this result for arbitrary delay distributions g{s). 
The stability of the solutions of Eq. Q is determined by 
linear stability analysis. We add a perturbation qk(t) to 




Global frequency of phase-locked steady states. We look 
for phase-locked synchronized solutions of Eq. with no 
phase lags: 
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FIG. 1: Examples of connection topologies considered in this 
work. Left to right: all-to-all connected system, ring topology 
with unidirectional connections, square lattice, and randomly 
connected system. Each oscillator has the same number of 
connections to other oscillators in all these systems. 
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h(t) = nt + eq k (i), 
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with e <§; 1 and substitute it into Eq. ([lj. The first 
order dynamic equations for the time evolution of the 
perturbation read: 



where 
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For the functional form of the perturbation we introduce 
the ansatz qk{t) = Cke xt and substitute into Eq. Mfy: 
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where g(X) is the Laplace transform of the delay distri- 
bution [TO [T2]: 
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Since dki = n from the requirement that any oscil- 

lator k has exactly n coupling neighbors, then for a =/= 
Eq. (TsT) can be rearranged: 



7(A) V" 



A 



N 



- + i)ck = J2 d 



kl Cl- 



(10) 



1=1 



Values of A such that g(X) = lead to singularities that 
are not meaningful physically, and must be excluded. We 
can write for the left hand side: 

n ( A \ , 
9(A) V" / 

where £ are the eigenvalues of the connectivity matrix D. 



In matrix notation Eq. (10) reads (c = Dc 



Using Gerschgorin's circle theorem [TPI [T3HT5] . an up- 
per bound on the modulus of the eigenvalues £ is given 
by the number n of connections each oscillator has: 
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from which: 
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The absolute value of the Laplace transform, |<?(A)|, is 
also bounded. With A = x + iy, and \g(s)\ = g(s), since 



g(s) > V s, we make use of the triangle inequality to 
show that for all x > 0: 



lff(A)| = 



dsg(s) e - [x+lv)s 
dsg{s)\e- xs \ \e tys \ < 1, 



< / ds g(s) 
o 



-{x J riy)s 



(14) 



where we have used that 
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[0, 1] V x > 0. Since g(X) and £ are complex num- 
bers, we define: 



g(X) = \g(X)\e^, 
C = IC|e l «, 



(15) 
(16) 



where ip and ^ are phases in the complex plane. Rewrit- 



ing Eq. ( 11 ) using this notation, we obtain: 

(x + iy + a), 
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which can be separated into real and imaginary parts: 

a\g(X)\^cos(ip + £) = x + a, (18) 
n 

a|KA)|^shH> + £) = y. (19) 
n 

Squaring these two last equations and adding them gives: 
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WithEqs. ( |13|14|18|20[ ) we can show the following result 
concerning the stability of the synchronized solutions of 
Eq. ([lj with frequency obeying Eq. Q: for all A that 
satisfy Eq. (10), Rc(A) < holds if and only if a = 
Kh'{—Q,f) > 0. This means that the stability of the 
synchronized states given by the solutions of Eq. Q 
depend only on the sign of a, which is independent of the 
shape of the delay distribution. This is the key finding 
of this paper. 

This result can be shown as follows. First we assume 
there exists an a > such that x = Re(A) > and show 
that this leads to a contradiction. We start dividing both 
sides of Eq. ((20l by a 2 : 
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Using from Eq. Q that |£| /n < 1, and from Eq. ( fli] ) 
that |ff (A) | < 1 V x > 0, we see that the product on the 
left hand side of Eq. ( 21 ) is smaller or equal than 1. The 



case of equality to 1 implies x = y = 0, i.e. A = 0. This 
corresponds to the eigenvector (1, 1, 1) T , which 
reflects the rotational symmetry of the system. The 
system is neutrally stable under such a perturbation. 
For all A 7^ 0, the right hand side of the equation will 
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FIG. 2: Left (dashed red) and right (continuous blue) hand 
side of Eq. (241 plotted as a function of x. The dashed curve 
always crosses the the j/-axis at y = — \a\ and the continu- 
ous curve always crosses in the interval [— \a\ , 0). Since the 
continuous curve approaches zero for x — > oo, but crosses the 
y-axis at values > — \a\, there must always be an intersection 
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FIG. 3: Re (A) versus flf for the eigenvalue A with the largest 
real part. Dashed blue line represents the discrete delay case 
and the solid orange line a distributed delay with an expo- 
of the two curves for x > 0, and hence a solution to Eq. (24}. nential distribution of mean f, g(s) = e~ s/T /f. N = 2, in- 
trinsic frequency ui = 0.223 min -1 , coupling strength K = 
0.07 min -1 , nearest neighbor interaction, and coupling func- 
tion h(9) = sin(6»). 

always be greater than 1. This poses a contradiction 
and it follows that for a > there are no solutions 
with Re(A) > 0. This means that a > assures linear 
stability. 

Now we complete the proof showing that if a < there 
is always at least one x = Re(A) > 0. For a < 0, we can 
rewrite Eq. (18) using a = — \ a\: 



and the right hand side of Eq. (24 1 



vertical axis at a value smaller than — \a\. 
for a < there always exists a solution x 
see Fig. [2j and therefore the system is linearly unstable 



cannot cross the 
It follows that 
Re(A) > 0, 
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A case differentiation for the cosine term is necessary to 
cover all possible signs of the left hand side. For the 
case cos(V> + £) < we have cos(ip + £,) = — \cos(ip + £)|. 
Rearranging Eq. ([22]) to find x yields: 



1 + |.g(A)| M |cos(^ + 01 ), (23) 



and we see that x > 0. For the case cos(V> + £ ) > 0, the 
equation 
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can be graphically analyzed plotting separately the left 
and the right hand sides, see Fig. [2] In the limit of x 
approaching infinity we have: x — > +oo |<?(A)| —¥ 
0; whereas for the case of x approaching minus infinity: 
x — > — oo => |.g(A) | — > 00, see Eq. (14). To finish the 



qualitative characterization of the right hand side of Eq. 
(24), we look at its behavior at x = 0. Recalling the 



normalization of the delay distribution and using that 



-iys 



1, we find: 



IS(A)| = 



ds g{s) e' iyt 



< 



/ dsg(s) |e-*"'| = 1, 
Jo 

(25) 



Conclusions and outlook. We have shown for systems 
of phase oscillators coupled with distributed delays that 
the linear stability of the phase-locked steady state with- 
out phase lags, given by the solutions of Eq. Q , depends 
only on the mean of the delay distribution, and is inde- 
pendent of its shape. A natural conclusion is that the 
discrete delay approximation is a reasonable one when 
interested in synchronized steady state dynamics, which 
have been the focus of much research. Among the many 
examples of studies using discrete delays, we highlight 
their use to model cellular coupling in the vertebrate seg- 
mentation clock [5]. This permitted the characterization 
of several mutations in zebrafish leading to the discov- 
ery of the first period mutants in this system [B]. Our 
work here shows that the results for the segmentation 
clock should be robust to the use of distributed delays, 
validating the discrete delay approximation in [BJ • By 
extension, it also validates this use in other problems [1] . 

Our result in this paper concerns steady states, but we 
can also say something about transient dynamics. Eq. 
(11) suggests that for different g(X) the real and imagi- 
nary parts of A might change, while the sign of the real 
part is independent of g(X). In this case, the modes re- 
sulting from perturbations of the steady state can decay 
or grow with different rates for delay distributions with 
different variance and skewness. This is indeed the case, 
as we show in Fig. |3j where we see that the fastest de- 
caying mode is different for discrete delay and a gamma 



distributed delay. Note that, despite the different shape 
of the curves for the different kind of delays, the crossings 
of the Re (A) = axis occur at the same points for both 
kind of delays, as our result assures. This dependence of 
the dynamics on the shape of the delay distribution is an 
intriguing topic that we will address in future work. 
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